
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
;
; The WRF ARW input file.  

  if (.not. isvar("input_file")) then
      start_date = systemfunc("date +%Y%m%d") + "00"
      input_file = systemfunc("ls /public/users/wrf/wrfout/" + start_date + "/wrfout_d01_*")
  end if

  if (.not. isvar("output_prefix")) then
      output_prefix = "/public/users/wrf/wrfpost/test/surface/surface_d01_" + start_date + "_"
  end if

  a = addfile(input_file + ".nc", "r")

  wks_type = "png"
  wks_type@wkWidth = 800
  wks_type@wkHeight = 800

; Set some basic resources
  res = True
  res@MainTitle                   = "WRF SURFACE"
  res@Footer                      = False
  res@NoHeaderFooter              = True

  pltres = True
  pltres@NoHeaderFooter       = True
  mpres = True
  mpres@NoHeaderFooter       = True
  mpres@mpDataBaseVersion       = "MediumRes"
  mpres@mpDataSetName           = "Earth..4"
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpGridLineColor             = "Black"
  mpres@mpLimbLineColor             = "Black"
  mpres@mpPerimLineColor            = "Black"
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpGridLineThicknessF        = 1.0
  mpres@mpLimbLineThicknessF        = 1.0
  mpres@mpNationalLineThicknessF    = 2.0
  mpres@mpUSStateLineThicknessF     = 1.0
  mpres@mpOutlineBoundarySets       = (/"National","Continental", "Geophysical"/)   ; turn on country boundaries
  mpres@mpOutlineSpecifiers = "China:states"
  mpres@mpGeophysicalLineColor      = "Black"       ; color of cont. outlines
  mpres@mpGeophysicalLineThicknessF = 2          ; thickness of outlines
  mpres@mpNationalLineThicknessF = 2          ; thickness of outlines




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Reading times
    char_times = a->Times
    hours = WRF_Times2Udunits_c(char_times, 0)
    ntimes = dimsizes(hours)         ; number of times in the file
    wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
    ; convert to Beijing Time
    hours = hours + 8
    times_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
    timesTideStr_cst = ut_string(hours, "%Y%N%D%H%M%S")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    contour_dt = 1 ; time interval of contour
    dt = tointeger(contour_dt / wrf_dt)
    do it = 0, ntimes-1, dt             ; TIME LOOP
    
    output_file = output_prefix + timesTideStr_cst(it)
    wks = gsn_open_wks(wks_type, output_file)

    print("Working on time: " + times_cst(it) )
    res@TimeLabel = times_cst(it)   ; Set Valid time to use on plots


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        

    slp = wrf_user_getvar(a,"slp",it)    ; slp
      wrf_smooth_2d( slp, 3 )            ; smooth slp
    tc = wrf_user_getvar(a,"tc",it)      ; 3D tc 
    td = wrf_user_getvar(a,"td",it)      ; 3D td 
    u  = wrf_user_getvar(a,"ua",it)      ; 3D U at mass points
    v  = wrf_user_getvar(a,"va",it)      ; 3D V at mass points
    td2 =  wrf_user_getvar(a,"td2",it)   ; Td2 in C
    tc2 = wrf_user_getvar(a,"T2",it)     ; T2 in Kelvin
      wrf_smooth_2d( tc2, 3 )
      tc2 = tc2-273.16                  ; T2 in C
      tc2@units = "C"
    u10 = wrf_user_getvar(a,"U10",it)    ; u at 10 m, mass point
    v10 = wrf_user_getvar(a,"V10",it)    ; v at 10 m, mass point

    tf2 = 1.8*tc2+32.                    ; Turn temperature into Fahrenheit
      tf2@description = "Surface Temperature"
      tf2@units = "F"
    td_f = 1.8*td2+32.                   ; Turn temperature into Fahrenheit
      td_f@description = "Surface Dew Point Temp" 
      td_f@units = "F"
;    u10 = u10*1.94386                    ; Turn wind into knots
;    v10 = v10*1.94386
;      u10@units = "kts"
;      v10@units = "kts"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ; Plotting options for T                
      opts = res                         
      opts@NoHeaderFooter       = True
      opts@cnFillOn = True  
      opts@cnInfoLabelOn     = False
      opts@ContourParameters = (/ -20., 40., 2./)
      opts@gsnSpreadColors = True
      opts@gsnSpreadColorEnd = -5  ; End third from the last color in color map
      opts@gsnSpreadColorStart = 5  ; End third from the last color in color map
      contour_tc = wrf_contour(a,wks,tc2,opts)
      delete(opts)


    ; Plotting options for Td
      opts = res         
      opts@cnFillOn = True 
      opts@cnLinesOn = True
      opts@cnLineLabelsOn = True
      opts@ContourParameters = (/ -20., 40., 2./) 
      opts@cnLineLabelBackgroundColor = -1
      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
;      contour_td = wrf_contour(a,wks,td_f,opts)
      delete(opts)


    ; Plotting options for SLP                     
      opts = res         
      opts@NoHeaderFooter       = True
      opts@cnInfoLabelOn        = False
      opts@cnLineColor = "Blue"
      opts@cnHighLabelsOn = True
      opts@cnLowLabelsOn = True
      opts@ContourParameters = (/ 900., 1100., 4. /)
      opts@cnLineLabelBackgroundColor = -1
      opts@gsnContourLineThicknessesScale = 1.5
      contour_psl = wrf_contour(a,wks,slp,opts)
      delete(opts)

    ; Plotting options for Wind Vectors                 
      opts = res         
      opts@NoHeaderFooter       = True
      opts@FieldTitle = "Wind"       ; overwrite Field Title
      opts@NumVectors = 47           ; density of wind barbs
      vector = wrf_vector(a,wks,u10,v10,opts)
      delete(opts)
  

    ; MAKE PLOTS                                       
      plot = wrf_map_overlays(a,wks,(/contour_tc,contour_psl,vector/),pltres,mpres)

    delete(wks)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

end
